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ABSTRACT 



This paper, and its companion, investigate the evolution of dense stellar 

Q ! systems due to the influence of two-body gravitational encounters, physical 

collisions and stellar evolution. Our goal is the simulation of the densest centers 

of galaxies, like M32, which reach stellar densities near 1O 6 M /pc 3 and which 

may harbor black holes. These systems have a different Safronov number (the 

dimensionless ratio of stellar binding energy to mean stellar kinetic energy) 

Q\ ! than globular clusters, substantially increasing the importance of physical 

collisions relative to gravitational encounters. In this paper, we focus only on 

the gravitational encounters. We demonstrate, first, that our simulations with 
00 

small N with a Hernquist tree code yield results basically in accord with years of 
effort studying globular clusters. Second, we investigate (crudely) core collapse 
in rotating systems with mass segregation, to separate out the effects purely due 
to two-body encounters from those seen in the more complex second paper. 
^ | The initial configuration for each experiment is an isotropic Kuzmin-Kutuzov 

model. Hernquist's tree code is utilized to simulate the dynamical evolution, 
with dynamical relaxation produced by the graininess of the low-N potential 
field. Consistent with previous studies, we find that systems whose constituent 
particles follow the Salpeter initial mass function rapidly undergo core collapse 
through Spitzer's mass segregation instability. All rotationally flattened 
systems show a decrease in flattening with time, consistent with Fokker-Planck 
calculations. 

An interesting new result concerns the well-established inability of simulators 
to identify a static center in their simulations of collisional systems. We find 
that the lump of high mass stars which condenses at the center wanders about 
the core in Brownian motion. We see this most clearly in the simulations with a 
mass function, in which the lump (possibly abetted by our force softening) lives 
on as a discrete entity after its formation. Even in the equal-mass case a similar 
but transient concentration appears and jitters about the center with similar 
radii to the unequal-mass case. 
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Subject headings: Galaxy: globular clusters: general - galaxies: nuclei - galaxies: 
star clusters 



1. Introduction 

There are many physical processes which govern the evolution of a dense stellar system. 
Two-body relaxation, stellar evolution and mass loss, star formation, and stellar collisions 
all influence the evolution of the system. We examine each of these processes through the 
use of N-body simulations. In order to discriminate among effects caused by these processes 
we have divided the study into two parts. This study (Paper I) examines the effects caused 
by two-body gravitational scattering, system rotation, and a stellar mass spectrum - the 
processes that do not depend on the finite size of the interacting particles. The purpose of 
this section of the work is twofold. First, we wish to demonstrate that our implementation 
of this technique yields results consistent with the extensive body of work on globular 
clusters, since the absence of such consistency would cast doubt on the more complicated 
results in paper II. Of course, we are also interested in any new effects that come out of 
this treatment of rotation, although they turn out to be rather meager. In Paper II we 
investigate the effects of stellar collisions, stellar evolution and mass loss, star formation, 
and a central black hole. Since several of the processes examined in Paper II either compete 
with or enhance mass segregation, it is necessary to examine them separately in order to 
determine their relative importance to the evolution of a dense stellar system. 

An examination of these processes may help to explain the structure of the dwarf 
elliptical M32, an E2 galaxy containing no dust or disk structure. Its proximity allows a 
detailed study of its kinematics. Hubble Space Telescope observations indicate a maximum 
core size of 0.37 pc and a central density in excess of 4 x 10 6 M pc -3 (Lauer et al. 1992). 
For comparison, the globular cluster M13 has corresponding values of r c = 1.8 pc and 
p c = 2.0 x 10 3 M pc" 3 (Djorgovsky 1993). The central velocity dispersion of M32 is 126 
km s _1 , and the rotation speed is 50 km s" 1 less than 4 pc from the center (van der Marel 
et al. 1997, 1994b). Dressier and Richstone (1988) have ruled out constant M/L models 
of M32 using a maximum entropy technique; we are left with two kinds of models: those 
which contain a central black hole, and those with a central cluster of either low-mass or 
degenerate stars (e.g. Goodman and Lee 1989, hereafter GL). Richstone et al. (1990) have 
shown that the observations are consistent with a central black hole of mass 0.7-3xl0 6 
M Q or a 10 7 M dark cluster with a core radius < 0."4 (1.3 pc). Subsequent studies (van 
der Marel et al. 1994a, Qian et al. 1995) showed that a good fit to the data is achieved with 
models containing a central black hole of mass of 1.8 x 10 6 M . Recently van der Marel et 
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al. (1997) have derived a mass of 3 x 10 6 M for the black hole. 

The time scales for two-body relaxation and physical stellar collisions place restrictions 
on the allowed core configurations of a dense system. The core dynamical time is given by 
(Binney and Tremaine 1987, hereafter BT) 



£dyn, c = \J3ir/l6Gp c , (1) 



and the core relaxation time is 



a 3 



t rehc = 0.337- -— - . (2) 

Gm*p c logA c 

Here p c is the core density, a is the one-dimensional velocity dispersion, m* is the stellar 
mass, and logA c is the Coulomb logarithm for the core. The time scale for a star in the 
core to suffer one collision is (BT) 



t, 



col,c 



16VW*(1 + 0)1 1 , (3) 



where n c is the stellar number density in the core and r* is the stellar radius. The Saffronov 
number is given by 

6 = ^- (4) 

Assuming a marginally unresolved 0.37 pc core of solar-mass stars, the central 
relaxation time is 6.5 x 10 7 y, and the collision time scale is 1.1 x 10 10 y. These time scales 
become 3.8 x 10 9 y and 6.5 x 10 11 y, respectively, 0.5 parsec from a 3 x 10 6 M black hole 
in the cusp model of Lauer et al. (1992). For comparison, the central relaxation time for 
M13 is 3 x 10 8 y and the collision time is 2 x 10 12 years. 

There are difficulties with both interpretations of the data. The short central relaxation 
time for the model lacking a central black hole implies that we have a privileged view to 
the core collapse of M32. If the collapse were reversed by the increased binding energy of 
binaries, or by the ejection of stars from the core, then M32 is in postcollapse reexpansion, 
and the "bounce" event cannot lie in the recent past (GL). The black hole interpretation 
suffers from what GL term the "luminosity problem" - the complete absence of any activity 
in the nucleus of M32. This is not necessarily fatal to the model either, if M32 is presently 
between episodes of flares caused by disrupted stars in the nucleus (Rees 1988). Presumably, 
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the only unique feature of M32 is its proximity, so that understanding its core structure will 
provide insight into the behavior of similar (but more distant) systems. We seek therefore 
to understand the processes which can lead to M32's remarkable nuclear properties. 

The evolution of stellar systems is usually studied through the use of Fokker-Planck 
(hereafter FP) calculations or N-body simulations. Each of these techniques possesses 
strengths and weaknesses. The principal advantage which FP calculations enjoy is their 
ability to handle large-N systems. Current N-body techniques can be employed to track 
up to 10 5 particles (e.g., Fukushige et al. 1991), an order of magnitude lower than globular 
clusters and several below small elliptical galaxies. This is a problem in that the various 
physical processes which dominate cluster evolution depend on different powers of N, so 
there is no simple way to scale the results from a small-N simulations to a large-N ensemble 
(e.g., Giersz and Heggie 1994). The small-N statistics are often too noisy for sound analysis. 
In addition, the dynamic range of FP calculations is excellent, whereas in many N-body 
codes it is limited by small N or the softening parameter approach to energy conservation. 

On the other hand, standard FP calculations make use of an orbit-averaged diffusion 
approximation to the system evolution. Thus discrete events, such as stellar collisions, 
evolution and mass-loss, must be treated statistically, whereas they are straightforward 
to incorporate into N-body codes. In addition, the distribution functions used in FP 
calculations are generally either functions of t and E, or at most t, E, and J, with no 
provision made for a possible third integral. N-body simulations do not suffer from this 
problem, since the distribution function need never be explicitly formulated in terms 
of conserved integrals. In addition, non-axisymmetric mass distributions are at present 
prohibitively complicated in an FP calculation, as are evolving off-center nuclear mass 
concentrations, whereas in an N-body simulation they require no extra labor. 

Nonetheless, there is quite good agreement between FP and N-body simulations for 
smaller-N systems under the influence of two-body relaxation. The two-body escape rate 
and core collapse occur in N-body simulations as predicted by FP calculations (Giersz 
and Heggie 1994), while provisions made to include large- angle scattering events in FP 
calculations have little or no effect on core collapse (Goodman 1985). These favorable 
comparisons make both techniques valuable for studying dense systems. 

An interesting aspect of evolution of a dense rotating stellar system is the secular 
evolution of its flattening. Although globular clusters are quite different from elliptical 
galaxies, the fact that the inner regions of some elliptical galaxies are dense enough to 
be relaxed makes globular clusters interesting laboratories wherein to study relaxation 
processes which may also operate in the cores of dense ellipticals. There is observational 
evidence to suggest that globular clusters become more spherical with age. Globulars 
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in the Magellanic Clouds, whose ages range from 10 7 to 10 10 y, may be somewhat flatter 
than Galactic globulars, which are thought to be 10 10 years old (Fall and Frenk 1985). It 
has been suggested that two-body gravitational relaxation in globular clusters transfers 
angular momentum to the halo, causing them to become more spherical. If true we expect 
this to hold for any rotating system whose relaxation time t r is less than the age of the 
universe T . This fits the general picture for globular clusters, where t r C T , and most 
elliptical galaxies, where t r 3> T . Agekian (1958) and Shapiro and Marchant (1976) used 
Maclaurin spheroids to model a rotating cluster with ellipticity e, and found that systems 
with e < 0.739 become more spherical as they evolve. Goodman (1983) showed that the 
Agekian model predicts that all flattened clusters become rounder with relaxation if there 
is a tidal truncation to the distribution function. He then analyzed cluster evolution using 
a single-mass FP calculation, and showed that cluster flattening does indeed decrease with 
time. Akiyama and Sugimoto (1989) ran 1000 equal-mass particles in an N-body simulation 
and showed that the angular momentum contained in cylindrical shells is transported away 
from the cluster center, with a characteristic timescale of 30 half-mass crossing times. 

As energy and angular momentum are transported out of the core it contracts and 
becomes denser, which in turn increases the transport rate, resulting in a shrinking 
isothermal core (Lynden-Bell and Wood 1968). Cohn (1980) followed this collapse over 20 
orders of magnitude in density using a single mass component FP calculation. Systems 
containing a range of masses evolve much more rapidly than systems with equal-mass stars 
(Spitzer and Saslaw 1966; Henon 1975; Chernoff and Weinberg 1990; Giersz and Heggie 
1996); thus we have examined the evolution of systems whose members have a range of 
masses as well as those containing stars of equal mass. 

One possible outcome of unchecked core collapse is an episode of stellar collisions 
resulting in the formation of a large black hole (Quinlan and Shapiro 1990) or a swarm 
of neutron stars, which in turn would collapse to form a black hole (Quinlan and Shapiro 
1989). Young (1980) computed models wherein a pre-existing black hole grows adiabatically 
at the center of a spherical cluster. Lee (1992) used FP calculations to study the evolution of 
rotating clusters containing massive black holes in their centers. We will delay a discussion 
of these topics until Paper II, where we treat stellar collisions, stellar evolution, and stellar 
mass loss, as well as the tidal disruption of low- J stars by a central black hole. 

It is often argued that the presence of binary stars will affect the evolution of a stellar 
system if the density of the system is sufficiently large (McMillan 1991 and references 
therein). The presence of a few "hard binaries" (i.e., binaries whose orbital speed is greater 
than the local velocity dispersion) in a globular cluster may be sufficient to halt core 
collapse (for a review see Spitzer 1987), but they are unequal to the task of supporting a 
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larger system like M32. A pair of 2M stars can supply an amount of energy equal to the 
total binding energy of a globular cluster. The binding energy of M32 is perhaps 5 orders 
of magnitude higher, however, and so binaries are not likely to affect the evolution of the 
core of M32. 



2. Initial simulation configuration 
2.1. Kuzmin-Kutuzov models 



We start each of our simulations with Kuzmin-Kutuzov (hereafter KK) models (Kuzmin 
and Kutuzov 1962; Dejonghe and De Zeeuw 1988). The models are either spherical, with 
zero rotation, or rotationally flattened to an axial ratio of 3.3:1. The potential-density pair 
for these models is 



<3>(g7, z) 



GM 

(zu 2 + z 2 + a 2 + c 2 + 2^d r c 2 + c 2 w 2 + ctz 2 ) 1 ' 2 



(5) 



and 



_ Mc 2 (a 2 + c 2 )w 2 + 2a 2 z 2 + 2a 2 c 2 + a 4 + 3a 2 Va 2 c 2 + c 2 w 2 + a 2 z 2 
P[W ' Z) ~ 4vr (a 2 c 2 + c 2 w 2 + a 2 z 2 f/ 2 (w 2 + z 2 + a 2 + c 2 + 2^a 2 c 2 + c 2 w 2 + a 2 z 2 f/ 2 ' ( ' 

Here zu, z, and suppressed are cylindrical radius, height, and azimuth. The ratio of the 
scale lengths, c/a, corresponds roughly to the axial ratio (see Dejonghe and De Zeeuw 1988 
for a discussion). 

The particle velocities are obtained with the KK density/potential pair and the Jeans 
equations. Assuming that the only non-zero first velocity moment is v^, and that the model 
is isotropic (i.e. a m = = a z ), cross terms like v m v z disappear from the Jeans equations, 
and we obtain 



_d_ 

dzu 



(pa 2 ) + p 




(7) 



and 



d (9$ 
d-z ipa) = ~ P dz- 



(8) 
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To evaluate a at the location of a specific particle we integrate the last equation over z 
(Satoh 1980; Binney et al. 1990), 



ro,z =oo 



pcr 2 (>,2:) = / P^jr-dz' (9) 

Jvj,z'=z UZ 

This can be used in to evaluate the rotation speed at a particle position in the first equation: 

vl 2 d , 9 . <9$ . . 

■uj ouj aw 

Figure [I] shows contours for v$ and a 2 for the KK model in which c/a = 0.3 and a + c = 1. 



2.2. Distribution sampling 

For a given value of c/a we use the density distribution given above and perform a 
rejection-based uniform sampling of a volume of space. Given a position (w, <fr, z), we put a 
particle at that location if 

p(w, z) . . 

, ! < R 11 
p(0,0) 1 ; 

where R is a random number on [0,1]. For each particle sampled at (zcr, z), we calculate 
the local dispersion and the particle's rotation velocity. The velocity assigned to the 
particle is simply \/3a, oriented randomly in space, plus v^4>. The resulting velocities in 
cylindrical coordinates for 3000 stars are shown in figure 0. In the spherical case the velocity 
distribution for all three projections is similar to the (—v m , v z ) plot. 

While this procedure assigns each star a speed on a velocity sphere, the stars phase mix 
in a very short time. It is important to note, however, that this procedure ensures that the 
simulated systems do in fact start out in global virial equilibrium with an isotropic velocity 
dispersion. As an example we show various kinematic quantities in cylindrical coordinates 
for the rotationally-flattened model II in Table [I]. Each of the averages is calculated 
globally. If we compare the binding energy V = —0.282551 to twice the total kinetic 
energy 2 x 0.129917 we get (2T + V)/\V\ = -8.04%. Although low-N noise contributes 
to the offset, the dominant source is the fact that a symmetric distribution about the 
mean nearest-neighbor distance produces an asymmetric binding energy distribution with a 
negative skew. 
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2.3. 



Initial mass function 



Two mass spectra were employed in these simulations. Let P(m) dm be the probability 
of a star's mass falling between m and m + dm. The first spectrum used was simply 
equal-mass stars. The other spectrum was a power law. 



For this case we used a = —2.35, the Salpeter initial mass function (Salpeter 1955). The 
spectrum is scaled by choosing a lower and an upper mass cutoff, m min and m max . 

The observational literature puts constraints on m min and m max . It is known that 
a main sequence star must have a mass of at least O.O8M in order for thermonuclear 
reactions to proceed in the core (e.g., Grossman, Hays, and Graboske 1974). The upper 
mass limit is somewhat more uncertain. The most massive star whose mass is directly 
measured (as a binary member) is HD47129; its mass is determined to be 3OM (Popper 
1980). An extrapolation of the mass- luminosity relation for an O star in 30 Doradus yields 
a mass of about 1OOM (Humphreys and Davidson 1986). 

For the simulations which use a Salpeter IMF we use m m i n = O.2M and 
mmax = 1OOM . Note that these are limits in the distribution function to be discretely 
sampled; an actual collection of 3000 objects will probably have an upper mass limit much 
less than 1OOM , since the probability density at the high-mass end is very low. If we have 
a mass spectrum characterized by N particles on an interval [m min , m max ], it is easy to show 
that the number of stars Nj on some smaller interval [mi, is 



We thus set m min such that we have at least a star or two on the subinterval [25, 100]M o 
for N=3000. Using m m j n =0.2 we expect Nj = 3.7 ± 1.9 stars on this interval. A typical 



The mass scale is arbitrary for the simulations in this study, since these systems evolve 
under the effects of two-body relaxation only. Although results are usually reported in units 
such that the total system mass is unity, we occasionally use solar masses to make explicit 
the observational constraints on the mass spectrum. This will also facilitate comparisons 
with the simulations of Paper II, wherein the various physical processes treated require 
knowledge of the mass scale. 



P(m) dm = A m a dm 



(12) 




(13) 



sampling of this distribution is shown in figure [| The maximum mass in this set is 32M . 
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Since stellar evolution and collisions were not allowed, the mass of each system and the 
number of stars were constant. Table |2] summarizes the model parameters for each of the 
simulations. In the second column a and c are the relevant scale lengths in the KK model. 
Column 3 refers to the distribution of stellar masses - either a collection of equal-mass 
stars or the Salpeter IMF - and columns 4 and 5 give the minimum and maximum masses, 
respectively, in the discrete sampling. Column 6 gives the total mass of the system in solar 
masses. 

2.4. Simulations 

To simulate the dynamical evolution of these systems we used the N-body tree code 
of Hernquist (1987, 1990), a Fortran implementation of the Barnes and Hut (1986) 
hierarchical force algorithm. The simulations were all performed on a Sun workstation. 
The output of the program consists of the mass, position, velocity, energy, and z angular 
momentum of each particle, written at specified intervals. All simulations conserved energy 
to better than 0.1%. 

3. Results 

The evolution in time of quantities presented in this section are plotted in units of 
the initial half-mass dynamical time. As it is instructive to gauge the progress of various 
processes with the core and half-mass relaxation time scales, these quantities are plotted in 
units of the half-mass dynamical time in figure f|. 

3.1. Mass segregation 

As Giersz and Heggie (1996) have noted, the most striking result of multi-mass 
model evolution is the rapid core collapse. Initially the heavy stars follow the same 
velocity distribution as the light stars, and so their temperature is much higher than their 
surroundings. They quickly give up energy to lighter stars and fall into the center of the 
system. It appears that this mechanism operates on roughly the half-mass relaxation time 
scale. In figure |5| we show the central region of model IV in z projection at four different 
times, in center-of-mass coordinates. A lump of high-mass stars quickly condenses near the 
center, ejecting low-mass stars. 

The rapid migration of heavy stars into the core appears to be the result of the mass 
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segregation instability (Spitzer 1969). For a system composed of two mass species (masses 
mi and m 2 , with mi < m 2 ; total mass in each species Mi and M 2 ), Spitzer (1969, 1987) 
defines the stability parameter x as 



for m 2 ^> rrii arid M 2 <C M\. If % < Xcrit = 0.16 then the system is stable against runaway 
mass segregation. If not, the self-gravity of the heavy stars requires a velocity dispersion 
that prevents equipartition with the lighter stars. The heavy stars give up energy to the 
lighter stars and sink to the center, while the light stars carry this energy into the halo. 
Since our simulations employed stars of many different masses we chose a separation mass 
m sep to categorize stars as either m x or m 2 . That is, if m* < m sep then we treat it as type 
mi, and if m± > m sep then we treat it as type m 2 . For every value of m sep chosen, x ^ XaAt- 
Thus our systems were unstable to runaway mass segregation. 

Following the work of King (1966) and the argument of Giersz and Heggie (1996) we 
define the "core radius" of the system as 



where a c and p c are the central velocity dispersion and density. Giersz and Heggie calculated 
the central quantities over the innermost 1% of the mass of the system, without worrying 
about the noise this incurs because of the simulation-averaging technique they employed. 
Since we ran only one simulation for each set of parameters we used the innermost 5% 
of the system mass in our core calculations. Although this produces slightly larger core 
radii, it is the behavior of the core radius defined at a particular level, and those quantities 
dependent upon it, which are of interest. 

Defined in this way, the core radius plummets in about a half-mass relaxation time, or 
10 half-mass crossing times, for multi-mass models III and IV (see figure ^). Once the core 
collapse phase has ended, the resulting core size depends on the softening length. Figure [7] 
shows the evolution of the core radius for softening lengths of 0.1 and 0.025. 

Our simulations show that core collapse slows as low-mass stars are evacuated from 
the core, consistent with the results of Giersz and Heggie (1996). Their depletion removes 
the dominant mechanism whereby energy is transported from the core to the halo (Spitzer 
1969, 1987). Figure |8| shows the z projection of model IV for r < 0.20 and r > 0.20, at 
t/tdjn,h = 102. Here r is measured from the lump center, defined as the mass center of the 




(14) 
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10 stars with the lowest potential energy. Global mass segregation is shown in figure |9], 
which plots the characteristic scale length (l/r) -1 for different mass groups as a function of 
time. 

Once the mass concentration forms it moves slowly around the central part of the 
system as defined by the stars not bound to it. The motion of the "density center" of a 
stellar system has been noted by several groups (although they have defined it in different 
ways). Miller and Smith (1992) described the behavior of the density center of their 
simulations as an "oscillation" and a "growing wave." Sweatman (1993), in contrast, argued 
that the motions are mainly a \/\fN noise effect. Spurzem and Aarseth (1996) Fourier 
analyzed the motion of the lump in their simulations and showed that it orbital timescale 
was about 14 half-mass crossing times (quite consistent with the timescale we observed). 
All of these studies involved systems of equal-mass stars. While we observe a similar 
phenomenon in our equal-mass simulations, it is much more striking in the mass-spectrum 
simulations. In these experiments the membership of the lump is nearly constant. 

A density enhancement also forms at the center of our equal-mass systems, but it has 
a somewhat more evanescent character, losing and accreting members. We believe that this 
is the same phenomenon, but it is certainly much less striking. It therefore seems likely 
that mass segregation plays an important role in the strength of the phenomenon. We 
believe that, whether long-lived or not, the behavior of the lump is Brownian motion due 
to scattering of individual stars. 



Figure [L(J shows the position of the model IV lump at 16 equally-spaced time intervals 
spanning the simulation. Using our r < 0.2 lump membership definition, we replaced all 
the lump members by a single particle with the former lump's mass and velocity. We 
then advanced the simulation by another 50 half-mass crossing times. The position of the 
lump replacement is shown in figure O. The initial and final positions, (0.11,0.10) and 



(0.19,0.61), respectively, have been circled. The initial velocity vector of the particle is in 
the positive x and negative y direction. 

The kinetic energy of the lump is 



-'lump g mm pl lumpl ■ \ 

We compare this quantity to the average kinetic energy per particle of the lump environment, 
defined as a cumulative average of the kinetic energies of particles not found within in the 
lump sphere: 
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-^V env 

E 



rriivf 



Tenv(r) = ^- . (17) 



A typical environment temperature profile is shown in figure [12]. T env is plotted as 
a function of distance from the system center-of-mass, for model IV at £/tdyn,h = 102. 
The lump temperature at this time is about 0.0011; model IV in general shows a lump 
temperature in the range 0.006 < T < 0.0012 and likewise for the environment temperature 
at low r. Although the low number of particles makes for noisy statistics, the evidence here 
suggests that the lump is in Brownian motion about the central region of the system, in 
thermal equilibrium with the immediate environment. 



3.2. Flattening 

An interesting aspect of the evolution of rotating systems is the change in the flatness 
of the system, as a function of both radius and time. One way to characterize this is to fit 
ellipses to projected light or density contours. The ellipticity is then defined as 



e=l-b/a, (18) 

where a and b are the major and minor semiaxes of the ellipse. Another way is to use 
the "dynamical ellipticity" (Goodman 1983), which characterizes the flattening as being 
due to rotation. It assumes that the three-dimensional density contours are concentric 
spheroids, and that the distribution function is a function of E and J only. While this 
definition is useful for FP calculations in which / = f(E,J), it is not appropriate for 
N-body simulations, where no assumptions about the distribution function are made. 

We seek to characterize the flatness of each system in a way that is numerically simple 
and that makes only minimal assumptions about the mass distribution. Defining 9 as 
latitude, we adopt 



e = l-(|tan0|), (19) 

where (| tan#|) is the density- weighted average of | tan#|. The flatness of a thick spherical 
shell bounded by r 1 and r 2 is thus defined as 
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T2 l"x/2 r2n 

tan#| p(r, 6, <fi) r 2 cos 9 dr d6< 



1 - JriJ - /2J ° , . (20) 

rT2 fir/ 2 f2ir v ' 

p(r, 6, <f>) r 2 cos 9 dr d6 dip 



hi J -it/2 Jo 

We now assume z axisymmetry and convert to a sum over discrete particles: 

N B he ll 

£ m,|tan^| 

' 1 -^7-, • W 

i 

Here the summations are performed over all particles found within the spherical shell of 
interest. For oblate systems e lies on (0, 1], while for prolate systems it spans [— oo, 0); for a 
sphere e = 0. Tests of this statistic can be found in the appendix. 

Consistent with the FP results of Goodman (1983), we find that our flattened models 
become secularly rounder. We plot ellipticity versus time for Lagrangian quintiles of model 
II in figure [T3] and of model IV in figure [TJ. It is clear that the system becomes significantly 
less flat in the inner 50% of the system. As expected, the decrease in flattening is more 
pronounced in the inner quintiles, where the local relaxation time scale is comparatively 
short. 

It appears that the reduction is flattening is reduced somewhat for the multi-mass 
model. While the flattening continues to decrease up to t = 280t dynjh in the equal-mass 
simulation, the reduction tails off around t = 150 — 200t dyn>h . This is probably due to 
the fact that light stars are ejected to the the halo rapidly during the early phases of the 
evolution as the heavy stars become centrally condensed; they gain energy more rapidly 
than they can lose angular momentum. 

A crude way to see the trend more clearly is to make a least-squares fit to the 



Lagrangian quintile ellipticities for all the models. Figure [15] shows this for the four models. 
In the flattened models, the outer 80% by mass of each system becomes rounder with time. 
The linear fits for the innermost quintile are not statistically significant for the Salpeter 
IMF models since the number of particles in the innermost quintile is low, resulting in 
excessive noise. 

It appears that the the inner region of the flattened equal-mass simulation becomes 
more spherical than the multi-mass simulation. We suspect this is due to a difference in 



the angular momentum transport in these regions mentioned above. In figure [16] we show 
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the angular momentum per unit mass for the inner three quintiles of the rotating models. 
Model II, the equal-mass system, clearly transports more angular momentum out of this 
region than does model IV. 

4. Summary 

In this study we have examined the evolution of dense stellar systems using N-body 
simulations. The results obtained concur with Giersz and Heggie (1996) - we find that 
systems whose stars have many different masses undergo rapid core collapse. The high-mass 
stars quickly condense into a dyamically distinct lump at the center. This core collapse 
terminates when the vast majority of low-mass stars, the carriers in the energy transport 
which produces core collapse, has been ejected from nucleus. The size of the lump is limited 
by the value of the softening parameter. We have determined that the lump, once formed, 
wanders about the the nucleus in Brownian motion, in thermal equilibrium with the stars 
in the nucleus of the system. 

We have devised a statistic to characterize the flatness of a system. We found that 
rotationally flattened systems become less flat in time, concurring with the FP results of 
Goodman (1983). This decrease in flattening is more pronounced in the inner regions of 
each system, where the relaxation time is comparatively short. We have found that a 
multi-mass system retains its flatness somewhat better than an equal-mass system over the 
inner half of their mass. This is probably related to the fact that the equal mass system 
transports more angular momentum out of its core than does the multi-mass system. 

This research was supported by NASA Theory Grant NAG 5-2758. 

A. Tests of the flattening statistic 

While the flattening statistic e defined in the text is computationally simple, it is not 
obvious that it is a good indicator of the "true" flatness of the system. Noise might be a 
problem in systems with only a few thousand particles, and the mathematical definition of 
e is formally singular if a particle crosses the z axis. In addition, while a radial modulation 
in density is necessary for any flattening to register, it will also introduce a systematic bias 
in the calculation of e. For example, we might expect our prescription to artificially inflate 
the value of the flattening if — d log p/d log r is large, since we would be biasing any mass 
sample with material at the equator. We investigated the validity of this prescription for 
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computing e in two ways. The first was an analytic calculation of e for an infinitesimally 
thin spherical shell in an ellipsoidal scale-free system. We assumed a density distribution of 
the form 



Po \^Jm 2 + a 2 z 2 J 
Here the parameter a is related to the intrinsic ellipticity ei nt by 



The calculated flatness is 



n/2 (cos 2 9 + a 2 sin 2 Of/ 2 



For even integral j3 > 2 we have 



(Al) 



6int = 1 - - • (A2) 

a 

The amount of material 5M contained in a thin shell of radius r and width 5r is 

2np r 3 5r W 2 cos 9d9 

SM = / , Z^TTKrTo • ( A3 ) 
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cos Odd 

_ ?r , 2 (cos 2 6»+a 2 sin 2 e)/ 3 / 2 



/9/2-1 

V i'(i-i)! _j 1 \ ntr ( a+V*^l \ 

Q 2 (2a)' 3 - 2 - 2 J(2i)! ^ a(2a)' 3 - 2 v / o^ : T 1U & U-v^^U 

1 - ; — • < A5 > 



Figure [I7| shows the calculated flatness versus the intrinsic ellipticity for three values of (3 
over the entire range of oblate objects. As expected, the bias introduced by this method is 
directed toward a greater flattening, although the amount is not very significant. In the 
outer regions of a KK model /3 — 4, so we do not expect significant errors. Even for the 
halo of a flattened Plummer law, where j3 = 5, the error would still not be intolerably large. 

Thus the second way we examined this statistic was through the use of Monte Carlo 
simulations. We sampled the power law mass distribution given by equation |A1| for 3000 
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particles outside of some minimum radius r min . We then divided the system into five 
equal-mass spherical shells and calculated the flattening of each. For each value of j3 and 
e int (setting r min = 1) we performed 100 experiments, finding the mean value of e and the 
standard deviation for each shell. The results are plotted in figure Consistent with 
our analytic calculations, the calculated flatness overestimates the intrinsic flattening, with 
the discrepancy increasing with increasing f3. Kuzmin-Kutuzov models decay as j3 = 4 or 
shallower, however, so the error is not large. 
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Table 1. Cylindrical velocity moments for the initial configuration of model II. This 
system is rotationally flattened to c/a = 0.3. 











coordinate 


Vi 






zu 


0.00317032 


0.04847563 


0.22014898 





0.32752717 


0.16131656 


0.23247047 


Z 


0.00576263 


0.05004266 


0.22362792 
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Table 2. 


Initial 


model parameters 


for the four N= 


=3000 simulations 


in this study. 


model 


1 — c/a 


N(m) 


"^smallest /M 


^largest /M 


M system /M 


I 


0.0 


S(m - 1M ) 


1.000 


1.000 


3000.0 


II 


0.7 


5(m - 1M ) 


1.000 


1.000 


3000.0 


III 


0.0 


m -2.35 


0.200 


32.01 


1890.7 


IV 


0.7 


m -2.35 


0.200 


32.01 


1890.7 
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g. 1. — Contours in and a 2 for a flattened Kuzmin-Kutuzov model. R plotted here is 
in the text. 
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Fig. 2. — Stellar velocities for 3000 stars sampled from a c/a = 0.3 Kuzmin-Kutuzov model. 
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Fig. 3. — A Salpeter spectrum sampled for 3000 bodies on [0.2,100]M G 




Fig. 4. — Core and half-mass relaxation times for the 4 models. All quantities are plotted 
in units of the initial half-mass dynamical time scale. 
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Fig. 5. — Z-projection of the core of model IV (t in units of the half- mass dynamical time). 




Fig. 6. — Core radius versus time for the four simulations. 
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Fig. 8. — z-projection of the core region of model IV for r < 0.2 and r > 0.2 at t/td yn ,h = 102. 
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Fig. 9. — Characteristic radii of different mass ranges for the multi-mass models. 



-30- 



(x,y) 



(y,z) 



(-X.Z) 



0.4 



I 1 1 1 I 1 1 1 
o 



o o 
o 

CP @ 
o o 



o 
o 



-0.4 



— I I I I I I I I I I I — 

-0.4 0.4 




o o 
o o 



O O ( 

-.0 o o 
3 Q> 

° o 



10. — Lump position at different times for model IV. 
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Fig. 11. — Position of model IV lump replacement particle. 
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Fig. 12. — Lump environment temperature for model IV. 




Fig. 13. — Flattening versus time for five Lagrangian quintiles of model II (equal-mass stars, 
c/a = 0.3). 




Fig. 14. — Flattening versus time for five Lagrangian quintiles of model IV (Salpeter IMF, 
c/a = 0.3). 




Fig. 15. — Linear least-squares fits to the ellipticities of the four simulations. 




Fig. 16. — Angular momentum per unit mass for the inner 3 Lagrangian quintiles of the 
rotating systems. 




Fig. 17. — Calculated thin-shell flatness versus intrinsic ellipticity for three different power 
law density profiles. 
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Fig. 18. — Calculated thick-shell flatness as a function of shell radius for (3 = 4, 5, 6. 
Within each f3 plot, the curves shown with one-a error bars are for, from top to bottom, 
e int = 0.8, 0.5, 0.2. The intrinsic flatness is plotted with dashed lines. 



